###############################################################################
#
#                            Abigail Hemmen, PhD Candidate
#                             University of Notre Dame 
#                                    08/29/2025
#                    Is it Worth Posting? Testing the Influence 
#                    of State Legislators in Simulated Social Media
#                          State Politics & Policy Quarterly
#
#                           R version 4.2.2 (2022-10-31 ucrt)
#                                 RStudio 2025.05.1 
#
###############################################################################
library(broom)
library(estimatr)
library(sandwich)
library(modelsummary)
library(MASS)
library(emmeans)
library(dplyr)
library(ggplot2)
library(sjPlot)
library(Hmisc)

# Set working directory for location of SPPQmain_clean.csv and control_RR_SPPQ.csv
setwd("C:/Users/12178/Dropbox/Dissertation/Why_Experiment/FINAL")

# Read in base data
data <- read.csv("SPPQmain_clean.csv")

data$treatment_simple <- factor(data$treatment_simple, levels = c("Control", "Policy", "Partisan", "Constituent"))

data$pid.1 <- factor(data$pid.1, 
                     levels = c("Democrats", "Independents", "Republicans"), 
                     ordered = F)

data$pid <- NA
data$pid[data$pid.1 == "Independents"] <- "I"
data$pid[data$pid.1 == "Democrats"] <- "D"
data$pid[data$pid.1 == "Republicans"] <- "R"
data$pid <- factor(data$pid,
                   levels = c("R", "I", "D"),
                   ordered = F)

########################## Models ####################################
### Approval Results
##Senator
approval_sen <- lm_robust(approval_sen ~ treatment_simple, 
                          data = data, se_type = "HC3")

#Copartisan
Coapproval_sen <- lm_robust(approval_sen ~ treatment_simple * copartisan, 
                            data = data, se_type = "HC3")
#Copartisan w/ covariates 
Controlsapproval_sen <- lm_robust(approval_sen ~ treatment_simple * copartisan + age + participation + sm_active,
                                  data = data, se_type = "HC3")

##Representative
approval_rep <- lm_robust(approval_rep ~ treatment_simple, 
                          data = data, se_type = "HC3")
#Copartisan
Coapproval_rep <- lm_robust(approval_rep ~ treatment_simple * copartisan, 
                            data = data, se_type = "HC3")
# Copartisan w/ covariates 
Controlsapproval_rep <- lm_robust(approval_rep ~ treatment_simple * copartisan + age + participation + sm_active,
                                  data = data, se_type = "HC3")


### Contact for potential outreach 
Fullopen <- lm_robust(open ~ treatment_simple, 
                      data = data, se_type = "HC3")
#Copartisans
Coopen <- lm_robust(open ~ treatment_simple * copartisan, 
                    data = data, se_type = "HC3")
#Copartisan w/ covariates
Controlsopen <- lm_robust(open ~ treatment_simple * copartisan + age + participation + sm_active,
                          data = data, se_type = "HC3")




### Trust
trust <- lm_robust(trust ~ treatment_simple, 
                   data = data, se_type = "HC3")
#Copartisans
Cotrust <- lm_robust(trust ~ treatment_simple * copartisan, 
                     data = data, se_type = "HC3")
#Copartisan w/ controls
Controlstrust <- lm_robust(trust ~ treatment_simple * copartisan + age + participation + sm_active,
                           data = data, se_type = "HC3")



################################## Exploratory Results ###############################

### Online activity 
SMsen <- lm_robust(approval_sen ~ treatment_simple * copartisan * sm_active, data = data, se_type = "HC3")

SMrep <- lm_robust(approval_rep ~ treatment_simple * copartisan * sm_active, data = data, se_type = "HC3")

SMopen <- lm_robust(open ~ treatment_simple * copartisan * sm_active, data = data, se_type = "HC3")

SMtrust <- lm_robust(trust ~ treatment_simple * copartisan * sm_active, data = data, se_type = "HC3")

### Political participation 
PARTsen <- lm_robust(approval_sen ~ treatment_simple * copartisan * participation, data = data, se_type = "HC3")

PARTrep <- lm_robust(approval_rep ~ treatment_simple * copartisan * participation, data = data, se_type = "HC3")

PARTopen <- lm_robust(open ~ treatment_simple * copartisan * participation, data = data, se_type = "HC3")

PARTtrust <- lm_robust(trust ~ treatment_simple * copartisan * participation, data = data, se_type = "HC3")

### Party Identification 
PIDapp_sen <- lm_robust(approval_sen ~ treatment_simple * copartisan * pid, 
                        data = data, se_type = "HC3")
PIDapp_rep <- lm_robust(approval_rep ~ treatment_simple * copartisan * pid, 
                        data = data, se_type = "HC3")
PIDopen <- lm_robust(open ~ treatment_simple * copartisan * pid, 
                     data = data, se_type = "HC3")
PIDtrust <- lm_robust(trust ~ treatment_simple * copartisan * pid, 
                      data = data, se_type = "HC3")


######################## Main Paper Text ##################################

########################### Main Effects ##############################

#### Approval 
#Senator
SenAppPolicy <- tidy(approval_sen, conf.int = TRUE) %>%
  filter(term == "treatment_simplePolicy") %>%
  mutate(model = "Policy Treatment", outcome = "Legislator 1")

SenAppPartisan <- tidy(approval_sen, conf.int = TRUE) %>%
  filter(term == "treatment_simplePartisan") %>%
  mutate(model = "Partisan Treatment", outcome = "Legislator 1")

SenAppConstituent <- tidy(approval_sen, conf.int = TRUE) %>%
  filter(term == "treatment_simpleConstituent") %>%
  mutate(model = "Constituent Treatment", outcome = "Legislator 1")

#Representative
RepAppPolicy <- tidy(approval_rep, conf.int = TRUE) %>%
  filter(term == "treatment_simplePolicy") %>%
  mutate(model = "Policy Treatment", outcome = "Legislator 2")

RepAppPartisan <- tidy(approval_rep, conf.int = TRUE) %>%
  filter(term == "treatment_simplePartisan") %>%
  mutate(model = "Partisan Treatment", outcome = "Legislator 2")

RepAppConstituent <- tidy(approval_rep, conf.int = TRUE) %>%
  filter(term == "treatment_simpleConstituent") %>%
  mutate(model = "Constituent Treatment", outcome = "Legislator 2")

# Combine all models
plot_data <- bind_rows(
  SenAppPolicy, SenAppPartisan, SenAppConstituent,
  RepAppPolicy, RepAppPartisan, RepAppConstituent
)

# Compute average estimates and average confidence intervals by treatment
avg_plot_data <- plot_data %>%
  group_by(model) %>%
  summarise(
    estimate = mean(estimate),
    conf.low = mean(conf.low),
    conf.high = mean(conf.high),
    .groups = "drop"
  )

# Ensure correct order for plotting
avg_plot_data$model <- factor(avg_plot_data$model,
                              levels = c("Policy Treatment", "Partisan Treatment", "Constituent Treatment"))

# Plot  ##### Figure 2 
ggplot(avg_plot_data, aes(x = model, y = estimate, color = model)) +
  coord_flip() +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  theme_minimal() +
  labs(
    y = "Average Effect on Legislator Approval",
    x = "",
    color = NULL
  ) +
  guides(color = "none") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))


# Create grid of predictor values
newdata <- expand.grid(
  treatment_simple = levels(data$treatment_simple),
  copartisan = c(0, 1)  
)

##### Approval - Representative
predictions_rep <- predict(Coapproval_rep, newdata = newdata, se.fit = TRUE)

pred_df_rep <- newdata %>%
  mutate(
    fit = predictions_rep$fit,
    se_fit = predictions_rep$se.fit,
    lower = fit - 1.96 * se_fit,
    upper = fit + 1.96 * se_fit,
    model_label = "Legislator 2"
  )

##### Approval - Senator
predictions_sen <- predict(Coapproval_sen, newdata = newdata, se.fit = TRUE)

pred_df_sen <- newdata %>%
  mutate(
    fit = predictions_sen$fit,
    se_fit = predictions_sen$se.fit,
    lower = fit - 1.96 * se_fit,
    upper = fit + 1.96 * se_fit,
    model_label = "Legislator 1"
  )

# Combine both datasets
pred_df_avg <- bind_rows(pred_df_rep, pred_df_sen) %>%
  group_by(treatment_simple, copartisan) %>%
  summarise(
    fit = mean(fit),
    se_fit = sqrt(mean(se_fit^2)),  # Pooled SE across the two models
    lower = fit - 1.96 * se_fit,
    upper = fit + 1.96 * se_fit,
    .groups = "drop"
  )



# Plot without legislator distinction
ggplot(pred_df_avg, aes(x = treatment_simple, y = fit,  #Figure 3
                        group = factor(copartisan),
                        color = factor(copartisan))) +
  coord_flip() +
  geom_point(position = position_dodge(width = 0.6), size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), 
                width = 0.2,
                position = position_dodge(width = 0.6)) +
  scale_color_manual(values = c("0" = "purple4", "1" = "khaki3"),
                     labels = c("Out-partisan", "Co-partisan")) +
  scale_y_continuous(limits = c(0, 0.9)) +
  labs(x = "Treatment Condition",
       y = "Average Legislator Approval",
       color = "Partisan Relationship") +
  theme_minimal(base_size = 14)


#### Trust 
predictions <- predict(Cotrust, newdata = newdata, se.fit = TRUE)

# Combine predictions with new_data
pred_df <- newdata %>%
  mutate(
    fit = predictions$fit,
    se_fit = predictions$se.fit,
    lower = fit - 1.96 * se_fit,  # 95% CI
    upper = fit + 1.96 * se_fit
  )

ggplot(pred_df, aes(x = treatment_simple, y = fit,  # Figure 4
                    group = factor(copartisan), color = factor(copartisan))) +
  coord_flip() +
  geom_point(position = position_dodge(width = 0.6), size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), 
                width = 0.2,
                position = position_dodge(width = 0.6)) +
  scale_color_manual(values = c("0" = "purple4", "1" = "khaki3"),
                     labels = c("Out-partisan", "Co-partisan")) +
  scale_y_continuous(limits = c(0, .4)) +
  labs(x = "Treatment Condition",
       y = "Average Trust in Government",
       color = "Partisan Relationship",
       title = NULL) +
  theme_minimal(base_size = 14)


# Table 2 uses models Trust, Cotrust, and Controlstrust comparing three separate models
# and using the same outcome variable
modelsummary(trust, output = "markdown", stars = T)
# Called "Simple" in the text
#treatment_simplePolicy = "Policy" in paper text
# treatment_simplePartisan = "Partisan"
# treatment_simpleConstituent = "Constituent" 
# N size = 1511 and p-values in markdown

modelsummary(Cotrust, output = "markdown", stars = T)# Called "Copartisan" in the text
#treatment_simplePolicy = "Policy"
# treatment_simplePartisan = "Partisan"
# treatment_simpleConstituent = "Constituent"
#copartisan = "Copartisan"
# treatment_simplePolicy:copartisan = "Policy x Copartisan"
# treatment_simplePartisan:copartisan = "Partisan x Copartisan"
# treatment_simpleConstituent:copartisan = "Constituent x Copartisan"
# N size = 1511 and p-values in markdown


modelsummary(Controlstrust, output = "markdown", stars = T)# Called Covariates in the text
#treatment_simplePolicy = "Policy"
# treatment_simplePartisan = "Partisan"
# treatment_simpleConstituent = "Constituent"
#copartisan = "Copartisan"
# treatment_simplePolicy:copartisan = "Policy x Copartisan"
# treatment_simplePartisan:copartisan = "Partisan x Copartisan"
# treatment_simpleConstituent:copartisan = "Constituent x Copartisan"
# age = "Age"
# participation = "Political Participation"
# sm_active = "Social Media Use Intensity"
# N size = 1511 and p-values in markdown

# Table 3 uses models Fullopen, Coopen, Controlsopen comparing three separate models
# and using the same outcome variable
modelsummary(Fullopen, output = "markdown", stars = T)# Called "Simple" in the text
#treatment_simplePolicy = "Policy" in paper text
# treatment_simplePartisan = "Partisan"
# treatment_simpleConstituent = "Constituent" 

modelsummary(Coopen, output = 'markdown', stars = T)# Called "Copartisan" in the text
#treatment_simplePolicy = "Policy"
# treatment_simplePartisan = "Partisan"
# treatment_simpleConstituent = "Constituent"
#copartisan = "Copartisan"
# treatment_simplePolicy:copartisan = "Policy x Copartisan"
# treatment_simplePartisan:copartisan = "Partisan x Copartisan"
# treatment_simpleConstituent:copartisan = "Constituent x Copartisan"
# N size = 1511 and p-values in markdown

modelsummary(Controlsopen, output = 'markdown', stars = T) # Called Covariates in the text
#treatment_simplePolicy = "Policy"
# treatment_simplePartisan = "Partisan"
# treatment_simpleConstituent = "Constituent"
#copartisan = "Copartisan"
# treatment_simplePolicy:copartisan = "Policy x Copartisan"
# treatment_simplePartisan:copartisan = "Partisan x Copartisan"
# treatment_simpleConstituent:copartisan = "Constituent x Copartisan"
# age = "Age"
# participation = "Political Participation"
# sm_active = "Social Media Use Intensity"



predictions <- predict(Coopen, newdata = newdata, se.fit = TRUE)

# Combine predictions with new_data
pred_df <- newdata %>%
  mutate(
    fit = predictions$fit,
    se_fit = predictions$se.fit,
    lower = fit - 1.96 * se_fit,  # 95% CI
    upper = fit + 1.96 * se_fit
  )

ggplot(pred_df, aes(x = treatment_simple, y = fit, # Figure 5
                    group = factor(copartisan), color = factor(copartisan))) +
  coord_flip() +
  geom_point(position = position_dodge(width = 0.6), size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), 
                width = 0.2,
                position = position_dodge(width = 0.6)) +
  scale_color_manual(values = c("0" = "purple4", "1" = "khaki3"),
                     labels = c("Out-partisan", "Co-partisan")) +
  scale_y_continuous(limits = c(0, .45)) +
  labs(x = "Treatment Condition",
       y = "Average Willingness to be Contacted",
       color = "Partisan Relationship",
       title = NULL) +
  theme_minimal(base_size = 14)




########################## Appendix ####################################

################# Appendix B 
#Table 6
print(pred_df_avg) # for appendix

#Table 7 uses approval_sen, Coapproval_sen, Controlsapproval_sen, 
#same labels as Tables 2 & 3
print(approval_sen)
print(Coapproval_sen)
print(Controlsapproval_sen)

#Table 8 uses approval_rep, Coapproval_rep, Controlsapproval_rep, 
#same labels as Tables 2 & 3
print(approval_rep)
print(Coapproval_rep)
print(Controlsapproval_rep)



SenAppPolicy <- tidy(approval_sen, conf.int = TRUE) %>%
  filter(term == "treatment_simplePolicy") %>%
  mutate(model = "Policy Treatment", outcome = "Senator")

SenPartisan <- tidy(approval_sen, conf.int = TRUE) %>%
  filter(term == "treatment_simplePartisan") %>%
  mutate(model = "Partisan Treatment", outcome = "Senator")

SenConstituent <- tidy(approval_sen, conf.int = TRUE) %>%
  filter(term == "treatment_simpleConstituent") %>%
  mutate(model = "Constituent Treatment", outcome = "Senator")

RepAppPolicy <- tidy(approval_rep, conf.int = TRUE) %>%
  filter(term == "treatment_simplePolicy") %>%
  mutate(model = "Policy Treatment", outcome = "Representative")

RepAppPartisan <- tidy(approval_rep, conf.int = TRUE) %>%
  filter(term == "treatment_simplePartisan") %>%
  mutate(model = "Partisan Treatment", outcome = "Representative")

RepAppConstituent <- tidy(approval_rep, conf.int = TRUE) %>%
  filter(term == "treatment_simpleConstituent") %>%
  mutate(model = "Constituent Treatment", outcome = "Representative")


# Combine the rows into one data frame
plot_data <- bind_rows(SenAppPolicy, SenConstituent, SenPartisan,
                       RepAppPolicy, RepAppConstituent, RepAppPartisan
) %>%
  mutate(
    outcome = factor(outcome, levels = c("Senator", "Representative")),
    model = factor(model, levels = c("Policy Treatment", "Partisan Treatment", "Constituent Treatment")))

# Figure 8 
ggplot(plot_data, aes(x = model, y = estimate, color = model, linetype = outcome)) +
  coord_flip() +
  geom_point(position = position_dodge(width = 0.4), size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                width = 0.2, position = position_dodge(width = 0.4)) +
  theme_minimal() +
  #geom_hline(yintercept=0) +
  labs(
    y = "Estimated Treatment Effect",
    x = NULL,
    color = "Rhetorical Type",
    linetype = "Account"
  ) +
  #ggtitle("Treatment Effects on Legislator Approval") +
  scale_linetype_manual(values = c("Senator" = "solid", "Representative" = "dashed")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

 
# Two different legislator approvals accounting for copartisan relationship

# Create grid of predictor values
newdata <- expand.grid(
  treatment_simple = levels(data$treatment_simple),
  copartisan = c(0, 1)  
)

##### Approval - Representative
predictions_rep <- predict(Coapproval_rep, newdata = newdata, se.fit = TRUE)

pred_df_rep <- newdata %>%
  mutate(
    fit = predictions_rep$fit,
    se_fit = predictions_rep$se.fit,
    lower = fit - 1.96 * se_fit,
    upper = fit + 1.96 * se_fit,
    model_label = "Legislator 2"
  )

##### Approval - Senator
predictions_sen <- predict(Coapproval_sen, newdata = newdata, se.fit = TRUE)

pred_df_sen <- newdata %>%
  mutate(
    fit = predictions_sen$fit,
    se_fit = predictions_sen$se.fit,
    lower = fit - 1.96 * se_fit,
    upper = fit + 1.96 * se_fit,
    model_label = "Legislator 1"
  )

# Combine both datasets
pred_df_combined <- bind_rows(pred_df_rep, pred_df_sen)

# Plot with shape to distinguish models
# Figure 9
ggplot(pred_df_combined, aes(x = treatment_simple, y = fit, 
                             group = interaction(copartisan, model_label),
                             color = factor(copartisan),
                             shape = model_label)) +
  coord_flip() +
  geom_point(position = position_dodge(width = 0.6), size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), 
                width = 0.2,
                position = position_dodge(width = 0.6)) +
  scale_color_manual(values = c("0" = "purple4", "1" = "khaki3"),
                     labels = c("Out-partisan", "Co-partisan")) +
  scale_shape_manual(values = c("Legislator 2" = 16, "Legislator 1" = 17)) + # circle and triangle
  scale_y_continuous(limits = c(0, .9)) +
  labs(x = "Treatment Condition",
       y = "Average Approval",
       color = "Partisan Relationship",
       shape = "Legislator") +
  theme_minimal(base_size = 14)

################# Appendix D
sum(data$mc1) #1268 #85%

sum(data$mc2) #942 #63%
################# Appendix E 

plot_model(PIDtrust, type = "pred", terms = c("treatment_simple", "pid", "copartisan"),
           title = "Trust in Government by Treatment, PID, and Copartisanship",
           axis.title = c("Treatment", "Trust")) #Nothing here
plot_model(PIDopen, type = "pred", terms = c("treatment_simple", "pid", "copartisan"),
           title = "Willingness to be Contacted by Treatment, PID, & Copartisanship",
           axis.title = c("Treatment", "Willingness to be Contacted")) #Nothing here
plot_model(PIDapp_sen, type = "pred", terms = c("treatment_simple", "pid", "copartisan"),
           title = "'Senator' Approval by Treatment, PID, & Copartisanship",
           axis.title = c("Treatment", "Approval")) #Almost something here
plot_model(PIDapp_rep, type = "pred", terms = c("treatment_simple", "pid", "copartisan"),
           title = "'Representative' Approval by Treatment, PID, & Copartisanship",
           axis.title = c("Treatment", "Approval"))


################# Appendix F
plot_model(PARTtrust, type = "pred", terms = c("treatment_simple", "participation", "copartisan"),
           title = "Trust in Government by Treatment, Participation, & Copartisanship",
           axis.title = c("Treatment", "Trust")) 
plot_model(PARTopen, type = "pred", terms = c("treatment_simple", "participation", "copartisan"),
           title = "Willingness to be Contacted by Treatment, Participation, & Copartisanship",
           axis.title = c("Treatment", "Willingness to be Contacted")) 
plot_model(PARTsen, type = "pred", terms = c("treatment_simple", "participation", "copartisan"),
           title = "'Senator' Approval by Treatment, Participation, & Copartisanship",
           axis.title = c("Treatment", "Approval")) 
plot_model(PARTrep, type = "pred", terms = c("treatment_simple", "participation", "copartisan"),
           title = "'Representative' Approval by Treatment, Participation, & Copartisanship",
           axis.title = c("Treatment", "Approval")) 

################# Appendix G
plot_model(SMtrust, type = "pred", terms = c("treatment_simple", "sm_active", "copartisan"),
           title = "Trust in Government by Treatment, Social Media Use, & Copartisanship",
           axis.title = c("Treatment", "Trust"),
           legend.title = "Social Media Use") 
plot_model(SMopen, type = "pred", terms = c("treatment_simple", "sm_active", "copartisan"),
           title = "Willingness to be Contacted by Treatment, Social Media use, & Copartisanship",
           axis.title = c("Treatment", "Willingness to be Contacted"),
           legend.title = "Social Media Use") 
plot_model(SMsen, type = "pred", terms = c("treatment_simple", "sm_active", "copartisan"),
           title = "'Senator' Approval by Treatment, Social Media Use, & Copartisanship",
           axis.title = c("Treatment", "Approval"),
           legend.title = "Social Media Use") 
plot_model(SMrep, type = "pred", terms = c("treatment_simple", "sm_active", "copartisan"),
           title = "'Representative' Approval by Treatment, Social Media Use, & Copartisanship",
           axis.title = c("Treatment", "Approval"),
           legend.title = "Social Media Use") 

################## Appendix I
data <- read.csv("SPPQmain_clean.csv")

old <- subset(data, select = c('state','pid.1','copartisan', "approval_sen",
                                 "approval_rep", "open", "trust", "treatment_pid",
                                 "treatment_simple")) #More variables than I'll need but just in case

old <- subset(old, treatment_simple == "Control") # Only comparing the different control options

old$approval1 <- old$approval_sen # This just makes the col names identical
old$approval2 <- old$approval_rep

new <- read.csv("control_RR_SPPQ.csv")
new <- subset(new, select = c('state','pid.1','copartisan', "approval1",
                              "approval2", "open", "trust", "treatment_pid",
                              "treatment_simple"))

old$version <- "Original"
new$version <- "New"

combined <- bind_rows(old, new)

#Cleaning the labels up 
combined <- combined %>%
  filter(treatment_simple %in% c("Control", "Environment", "True")) %>%
  mutate(control_type = case_when(
    version == "Original" & treatment_simple == "Control" ~ "Original Control",
    version == "New" & treatment_simple == "Environment" ~ "Social Media Control",
    version == "New" & treatment_simple == "True" ~ "True Control"
  ))


####### Means Comparisons / Marginal Means 
ggplot(combined, aes(x = control_type, y = approval1)) + # Figure 22
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.2) +
  labs(
    x = "Control Condition Type",
    y = "Approval",
    title = "Comparison of Mean Approval of Legislator 1 by Control Group Version"
  ) +
  theme_minimal()

ggplot(combined, aes(x = control_type, y = approval2)) + # Figure 23
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.2) +
  labs(
    x = "Control Condition Type",
    y = "Approval",
    title = "Comparison of Mean Approval of Legislator 2 by Control Group Version"
  ) +
  theme_minimal()


ggplot(combined, aes(x = control_type, y = trust)) + # Figure 25
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.2) +
  labs(
    x = "Control Condition Type",
    y = "Trust",
    title = "Comparison of Mean Trust by Control Group Version"
  ) +
  theme_minimal()

ggplot(combined, aes(x = control_type, y = open)) + # Figure 24
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.2) +
  labs(
    x = "Control Condition Type",
    y = "Willingness to be Contacted",
    title = "Comparison of Mean Willingness to be Contacted\n by Control Group Version"
  ) +
  theme_minimal()



